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Monte Carlo methods play important part in modern statistical physics. The applica- 
tion of these methods suffer from two main difficulties. The first is caused by the relatively 
small number of particles that can participate in any numerical calculation. This means 
that scales larger than or comparable to the one that can be simulated are not taken into 
account. The second difficulty is the locality of the conventional Monte Carlo algorithms 
which leads to very ( sometimes unreasonably) long equilibration times. These obstacles 
can be eliminated in the framework of the multilevel Monte Carlo method described here. 
The basic approach is to describe the system at increasingly coarser levels defined on 
increasingly large domains, and transfer information back and forth between the levels 
in order to obtain a selfconsistent result. The method is illustrated for a test case of 
one-dimensional fluids 

1. INTRODUCTION 

The modern theory of classical liquids is based on the statement that the macroscopic 
characteristics of a many-particle system can be obtained by averaging over microscopic 
configurations, with the probability proportional to the Gibbs distribution function. Every 
microscopic configuration is defined by assigning specific locations to the particles, and the 
main computational problem of statistical physics consists in the difficulty of averaging 
over the enormous space of possible configurations. In order to estimate the value of this 
average the Monte Carlo technique for the canonical ensemble was proposed pp. 

However, the straightforward application of this technique is restricted to a small volume 
of the system under consideration, because one can consider only a relatively small number 
of particles in any numerical calculation. In order to minimize the surface effect on the 
bulk values which are to be calculated, artificial periodic boundary conditions are supposed 
PP similar to the case of an infinite crystal [2]. 

This means that at scales comparable with the periodicity cell size the fluctuations 
of the particle number are cut off. In order to avoid this difficulty several approaches 
have been suggested using the grand canonical ensemble The main trait of 

these algorithms is the generation of density fluctuations in the basic cell by adding and 
deleting particles in accordance with the value of the chemical potential. 

The Monte Carlo simulations both in canonical and grand canonical ensembles are 
very local, moving, e.g., one particle at a time. The operation of adding/deleting a 
particle is also local and needs to be followed by a local equilibration. This leads to 
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very slow changes of large-scale features, such as averages and various types of clusters 
(regions of (anti)aligned dipoles, crystallized segments, etc.). The larger the scale the slower 
the change and longer (per particle!) is the Monte Carlo process required to produce 
new independent features. Since many independent features are needed for calculating 
accurate averages, and since very-large-scale features need to be sampled, especially in 
the vicinity of phase transitions, the computations often become extremely expensive, 
sometimes even losing practical ergodicity. In recent years, a number of novel Monte 
Carlo algorithms ensuring ergodicity were proposed [HI, [Zl, |H]- In order to avoid the 
slowness of the ordinary Monte Carlo simulation, steps of more collective nature are used. 
Nevertheless, the simulation is performed in the periodicity cell and the estimation of 
bulk values may only be done by extrapolation p!^. 

An approach which allows to simultaneously overcome slowness and finite size effects of 
the conventional Monte Carlo method consists of a multilevel view of the system, realized 
by multilevel algorithms [HI, PUj. The multilevel algorithms construct a sequence of 
descriptions of the system under consideration at increasingly coarser levels and transfer 
information back and forth between the levels in order to obtain a selfconsistent result. 
The efficiency of multilevel methods in solving problems of statistical physics has been 
shown on examples with sufficiently simple systems jllj. It follows from these results that 
the effect of slowing down can be eliminated. Moreover, due to the possibility to simulate 
large volumes at coarse levels, the volume factor (the proportionality of the computational 
cost to the volume being simulated) can be suppressed as well. This means that the scale 
of treating a system is not restricted by the size of the basic cell. The elimination of both 
the slowing down and the volume factor allows one to investigate in the framework of the 
multilevel method long range phenomena such as phase transition. This possibility has 
been shown on two-dimensional examples of variable-coefficient Gaussian models jTTj and 
the Ising model P^. 

The successful application of the multilevel methods to lattice systems excites interest in 
adapting them to more complicated cases. The present paper treats multilevel algorithms 
to the investigation of fluids. Numerical results of simulations for one dimensional models 
are presented. These models reproduce common properties of fluids and allow comparisons 
to exact solutions. 

2. MULTILEVEL MONTE CARLO APPROACH 

2.1. Conventional Monte Carlo Method 

The Monte Carlo method in the statistical theory of liquids is used to evaluate numer- 
ically the average A of any functional A, defined by: 



where w{2Q is the probability density of the state X in the configuration space Q, and 
the nodes Xi are generated by a random walk in Vt that satisfies detailed balance. 




(1) 
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The simplest definition of tlie probability to walk from X to 2C_ in detailed balance is 
given by : 



uj{X ^20 = min 



(2) 



provided that the probability -P(X^?|X) for having chosen the state 2£L as the candidate 
to replace the current state X is symmetric, i.e., P(X'7\X) = P(X7\X'). 

The kind of probability density given by Gibbs in the canonical ensemble for simple 
liquids is: 

w{2Lj = const ■ exp( — ) (3) 

Kb ■ T 

where fi is the location of the i-th particle (i = 1, . . . , N), ks is the Boltzmann constant, 
T is the temperature and U is the potential energy. It is usually assumed to have the 
form: 

U = j:<l>i\rr-r,\) (4) 

i<j 

where 0(| fj — fj |) corresponds to the energy of a two-body interaction. 

In any numerical calculation it is possible to consider only relatively small number of 
particles (extremely small in comparison with Avogadro number). In order to minimize 
the surface effect in a small simulation volume, periodic boundary conditions are supposed 
pp. This means that the space is filled by particles which are located at points: 

R.^^=fi + f-L, l<i<N (5) 

where / is a vector whose components are integer, particles at points are located inside 
the basic cell of the linear size L, and N is the particle number in the periodicity cell. 

It follows from the periodicity condition that the real system is replaced by a super- 
lattice with the same value of the particle number density p = N/ (where D being the 
dimension of the space) in each cell. 

The motion of particles is continuous and in any configuration the particles have arbi- 
trary locations. The transition between states is made by the shift of one particle at a 
time P] by a small amount via the following displacement: 

fr'" = r,'''' + Sma.-{l-2-0 (6) 

where fj°''^ and r^"'^"' are old and new locations of the i-th particle, Smax is the maximum 
possible displacement of the particle along a coordinate axes, 1 is the vector all of whose 
components are 1, and ^ is a vector whose components are random numbers distributed 
uniformly on the interval [0, 1]. 

For shifting the particle with the number i, say, one can see from (0) that it is enough 
to use, instead of the Gibbs function with the energy (0}, the con(iztzona/ probability 
defined by: 

P{ri I Ri) = const ■ exp(— nj(Ri)/^5 T) (7) 
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The last relation gives us the probability to find the particle with the number i at the posi- 
tion when the locations of all other particles, defined by the set Ri = {ri, . . . , rj_i, rj+i, . . . , tn}, 
are fixed. The one-particle energy Mi(Ri) in ^ is defined by: 

«,(Ri)= ^ ^(Irl-r-;- I) (8) 

It follows from © that the conventional Monte Carlo simulation is the local process. 
Therefore one can expect the high speed convergence only in the case of systems with 
short-range correlations. 

2.2. Coarse variables 

A slowing-down is inherent not only in the conventional Monte Carlo algorithm, it is 
a common problem for all local processes (e.g. Gauss- Seidel relaxation for discretizied 
partial differential equations). The solution to this problem lies in introducing system 
changes of more collective nature. In the case of partial differential equations fast conver- 
gence of solutions had been attained by multigrid algorithms ^2)- These algorithms are 
looking for solutions on a sequence of lattices with increasingly larger meshsize (coarser 
scales) by combining local processing at each scale with various inter-scale (inter-lattice) 
interactions. A similar technique can be applied to the simulation of liquids. 

A possible way to introduce a coarse description of liquid consists in the discretization 
of space. The periodicity cell is divided into M disjoint parts (e.g. cubes) of equal 
volume with linear size hi, 1 < i < M (each being associated with a gridpoint i of 
the first coarse- level lattice). Configurations of the finest (particle) level are mapped to 
the first coarse level by the operation of coarsening, this operation creates the coarse-level 
variable set. For example, at any instant the corresponding coarse-level variable can be 
defined by coarsening the particle number: 

nl = Number of particles in V^^ (9) 

M 

with J2 = N, where is the total number of particles in the periodicity cell. The 

1=1 

set {n}} defines the current configuration on the first coarse-level: instead of particle 
locations the occupation numbers at gridpoints are used. 

Generally, the aim of coarsening is the creation of configurations represented by collec- 
tive variables which describe collective particle motions at different scales. The variable 
at each lattice point at each coarse level is defined as a local spatial average (an average or 
actually a sum over a certain neighborhood of the lattice point) of similar variables at the 
next finer level. The total value of each such variable is well defined for each configuration 
of particles given at the finest level. 

The extension of the coarsening operation © to coarser levels leads to the following 
definition of the coarse- variable at the level k: 

n) = k > 1 (10) 

for each volume element Vj' of level k, assuming it to be a union of volume elements of the 
level k — 1. The coarsening can be repeated till the coarsest level, whose choice depends 
on the scale of the phenomena one wants to compute. 
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The coarser level simulations are performed in periodicity cells of larger sizes. In this 
way finite size effects are suppressed and it is possible to reach the macroscopic description 
of the system. No slowdown should occur provided that the coarsening ratio (the ratio 
between a coarse meshsize and the next finer meshsize), as well as the average number of 
of original particles per mesh volume of the finest lattice, are suitably low. The typical 
meshsize ratio is 2, typical number of particles per finest lattice mesh is between 2 and 
10 (being usually larger at a higher dimension). More aggressive coarsening ratios would 
require much longer simulations to produce information for coarser levels. 

There are many possible ways to choose the set of coarse variables. A general criterion 
for the quality of this set is the speed of equilibration of a compatible Monte Carlo (CMC). 
By this we mean a Monte Carlo process on the fine level which is restricted to the subset 
of fine-level configurations whose local spatial averages coincide with a fixed coarse-level 
configuration. A fast CMC equilibration implies that up to local processing all equilibrium 
configurations are fully determined by their coarse-level representations (their local spatial 
averages), which is the main desired property of coarsening. 

The compatible Monte Carlo equilibration speed can be tested at the fine level, in the 
following way: after local thermalization of the ordinary Monte Carlo process, for a given 
configuration in the equilibrium the corresponding coarse variables are calculated by © or 
(jlUp . Then an initial configuration in the original variables is created in accordance with 
the coarse variables set. A following ensemble of compatible Monte Carlo processes yields 
an accurate estimate of the rate of approach to the equilibrium. The interpolation from 
the first coarse level configuration to the finest (particle) level needs a special consideration 
owing to the mapping of variables defined on a grid into the continuous space of particle 
coordinates. 



2.3. CP Tables and Multilevel Algorithm 

The main idea of the multilevel approach is to equilibrate on each level only modes 
with short (comparable with the meshsize) wave lengths. Long wave modes with slow 
convergence at a given level are equilibrated at coarser levels where their wave lengths 
are comparable with the meshsize. As a result, the multilevel process leads to fast equi- 
libration of all modes. 

In the framework of this multilevel Monte Carlo algorithm, only a local process is per- 
formed at each level, defined in terms of the corresponding variables. In order to calculate 
transition probabilities (j21) conditional probabilities similar to (|7j) should be derived for 
each coarse level. Such conditional probabilities are expressed in the form of a Conditional 
Probability (CP) table, which in principle tabulates numerically the probability distribu- 
tion of any pair of neighboring coarse-level variables given the values of all others. Of 
course, not a// other variables should in practice be taken into account: only the immedi- 
ate neighborhood of gridpoint pair under test counts, due to the near locality property of 
the conditional probability (cf. also discussion of near locality in [T^). 

For example, in terms of variables Q, (fTUI) defined at gridpoints a conditional probabil- 
ity tables Pk{n^, | sj, . . . , s-) can be constructed from a given sequence of configurations 
in equilibrium on the next finer level k — 1. These tables give us the dependence of the 
probability to find and particles at the two neighboring gridpoints i and j on / 
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values in their neighborhood by: 
q G neighborhood of i and j 



s..= Y: l<m<l (11) 



where are preassigned, suitably chosen coefficients, with a™' = a^" = for all m. 

For example, in one dimension j = i+1 and a possible choice, for / = 1, a}_i = a}_^_2 = 1, 
otherwise a™ = 0. 

In the coarse level Monte Carlo run, each trial move consists of particle exchange 
between two neighboring gridpoints, i.e. n^ — *■ nf = + An, Uj —>■ nj' = nj — An. In 
accordance with Q, the acceptance probability for this move is: 



ujiK^ X') = min 



Pki 


nf. 


hi 1 1 

nJ \sl. 




Pk 




k 1 1 

nJ 1 si, . . 





(12) 



The CP tables for any coarse level k are calculated by gathering appropriate statistics 
during Monte Carlo simulations at the next finer level k — 1. Because of the near-locality 
property, no global equilibration is needed; local equilibration is enough to provide the 
correct CP values for any neighborhood for which enough cases have appeared in the 
simulation. Thus, the fine-level simulation can be done in a relatively small periodicity 
cell. 

However, since the fine-level canonical ensemble simulations use only a small periodicity 
cell, many types of neighborhoods that would be typical at some parts of a large volume 
(e.g., typical at parts with average densities different than that used in the periodicity 
cell) will not show up or will be too rare to have sufficiently accurate statistics. Hence, 
simulations at some coarse level may run into a situation in which the CP table being used 
has flags indicating that values one wants to extract from it start to have poor accuracy. 
In such a situation, a temporary local return to finer levels should be made, to accumulate 
more statistics, relevant for the new local conditions. 

To return from a coarse level to the next finer level one needs ffist to interpolate, i.e., to 
produce the fine level configurations represented by the current coarse level configuration, 
with correct relative probabilities. The interpolation is performed by CMC sweeps at the 
fine level (few sweeps are enough, due to the fast CMC equilibration). This fast equili- 
bration also implies that the interpolation can be done just over a restricted subdomain, 
serving as a window: In the window interior good equilibrium is reached. Additional 
passes can then be made of ordinary (not compatible) MC, to accumulate in the interior 
of the window the desired additional CP statistics, while keeping the window boundary 
frozen (i.e., compatible). The window can then be coarsened (by the local spatial av- 
eraging) and returned to the coarse level, where simulations can now resume with the 
improved CP table. 

Iterating back and forth between increasingly coarser levels and window processing at 
finer levels whenever missing CP statistics is encountered, one can quickly converge the 
required CP tables at all levels of the system, with only relatively small computational 
domains employed at each level. The size of those domains needs only be several times 
larger than the size of the neighborhoods being used (with a truncation error that then 
decreases exponentially with that size). However, larger domains are better, since they 
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provide sampling of a richer set of neighborhoods (diminishing the need for returning later 
to accumulate more statistics), and since the total amount of work at each level depends 
anyway only on the desired amount of statistics, not on the size of the computational 
domain. 



3. NUMERICAL TESTS 

Thermodynamic properties of a one- dimensional fluid can be described exactly under 
the assumption of the nearest neighbor finite range interparticle interaction One 
can cite as an example of the short-range potential the truncated Lennard- Jones potential: 




for r < Tcut (^]^3^) 
for r > Tcut 



where a is related to the diameter of atoms, e is the measure of interparticle interaction 
and Tcut is the cut-off distance. 

At high temperatures the attractive part of the potential can be neglected and a good 
approximation for the interparticle interaction is the hard rods model: 

_MW r , r < 



1 .r>d 



where d is the diameter of every particle. 

In order to test the multilevel algorithm it was applied to the simulation of one- 
dimensional fluids. First of all, the speed of the compatible Monte Carlo equilibration was 
tested in the case of the interparticle interaction (fT!^ with cut = 2.5 ■ a. After locally equi- 
librating an ordinary Monte Carlo process, a sequence of equilibrium configurations was 
picked out. For each configuration the coarse- level variable set is defined by ® and an 
initial configuration for compatible Monte Carlo simulation is constructed in accordance 
with the frozen equilibrium coarse level structure. 

For the Lennard- Jones fluid a convenient quantity measuring the approach to thermal 
equilibrium is the energy of the system. However, this quantity is subject to fluctuations 
during a single run of the equilibration. In order to smooth out simulation results a number 
of compatible Monte Carlo simulations are performed using the same initial configuration 
but different sequences of random numbers. The relaxation function is defined by the 
averaging over these CMC runs: 

< Ujr) >-< UjTf) > 
= f/o- < U{rj) > ^''^ 

where U{t) is the energy of the system after sweep r, subscript / means final configura- 
tions, and Uq = U{0). 

The ensemble average relaxation function ()15|) is used to estimate the relaxation time: 
re,(f/o) = J Ruoir)-dT (16) 
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Here the relaxation time depends on the initial configuration of CMC runs, therefore addi- 
tional averaging of the relaxation time ()16|) is done over all the equilibrium configurations 
which were chosen during the ordinary Monte Carlo process as the initial configuration 
set for CMC: 

Teq =< TeqiUo) >Uo (17) 

Definitions and (|17p lead to the exact result in the case of a single classical relax- 
ation behavior. Nevertheless, even if the relaxation behavior is polydispersive, the integral 
(Uni) still gives an unambiguous measure of the relaxation rate [T7j. 




Figure 1. Calculated temperature dependence of the Compatible Monte Carlo energy 
relaxation time. The particle number density p-a = 0.2. o - 16 particles in the periodicity 
cell ( L = 80-cT, M = 32); A - 256 particles in the periodicity cell (L = 1280 -a, M = 512). 
The size of a subdomain was chosen to be /ii = 2.5 ■ cr. 



The temperature dependence of the relaxation time estimated in accordance with (jl7|) 
is shown in Fig^ In contrast to the ordinary Monte Carlo process the relaxation of CMC 
is not sensitive to the size of the periodicity domain and a fast equilibration implies that 
coarsening the particle number satisfies (in the case of the interaction potential (fT^ and, 
hence, in the high temperature limit ()14|) the requirements for the quality of the coarse 
variables set. 

In the special case of the one-dimensional fiuid, the main contribution to the distribu- 
tion of the particle number at any lattice point is expected to come from its nearest 
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neighbors. Therefore a simple possible form of the conditional probability table is given 
by: 



I ; '"J I ^jj) • • • ; ■^ijj — ^ kyil'i ; "'i+l I "'i+2) 

Quantities of were calculated at each lattice point i for each produced configuration at 
the next finer level, and from them the probability distributions (fT^ were accumulated, 
forming the CP table. Then, multilevel Monte Carlo simulation was performed using this 
CP tables. 

A suitable quantity for comparing simulation results at different levels is the fiuctuation 
z/fc of the particle number in the subdomain of size V'^ (corresponding to the gridpoint of 
the lattice with meshsize hk = hi ■ 2^^^): 



< > 



(19) 



where the mean values are defined by probabilities pkiji) to display the particle number 
in a volume of size V^: < >= ■ pk{n''), < >= ■ pk{n'^). 



> 0.4 




Figure 2. Particle number fiuctuation at different levels (hard rods fluid, h\ = 2.5 ■ d). 
Dashed lines correspond to exact bulk values (|2(jp. 



The fluctuation of the particle number does not depend on the volume if it is suffi- 
ciently large (as compared with the correlation length) and is related with the isothermal 

compressibility k = ( — ) UM- For ID hard core system the equation of state 

V / T= const 
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is found analytically (it is known as the Tonk's equation p!^) and the isothermal com- 
pressibility is given by: 

= p-keT ■K={1- p-df (20) 

Results of particle number fluctuation multilevel measurements in volumes of different 
sizes for hard rods systems are presented in FigI21 The dependence of this quantity on 
the subdomain size is explained by the size effects and coincides with the expression of 
the grand canonical ensemble |14j : 

C , , 

yk = yoo + ^ (21) 

where v^q is the bulk value of the particle number fluctuation ()20|) . The coefficient C as 
well as depend on the particle number density. One can see from FigOthat at high 
temperatures (the interaction model ()14|) five levels are enough in order to achieve the 
bulk behavior in subdomains connected with gridpoints of the coarsest level. 



pa=0.2 



5.5 




Figure 3. Particle number fluctuation at different levels (Lennard- Jones fluid, h\ = 2.5-0"). 



In the case of the Lennard- Jones fluid the particle number fluctuations depends also on 
the temperature. The dependence on the meshsize (the subdomain volume) for system 
with the interaction potential ()13p is shown in FiglS] for different temperatures. One 
can see that at high temperatures the properties of the Lennard- Jones fluid are similar 
to the hard rods system (see Fig|2I). At low temperatures the behavior of fluctuations 
drastically changes. The results obtained indicate that in contrast to the hard rods system 
the Lennard- Jones system loss homogeneity at low temperatures on fine scales. 
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4. CONCLUSION 

It follows from the present results that in the framework of the multilevel Monte Carlo 
method a suitable choice of coarsening allows to attain the macroscopic behavior. The 
particle number fluctuation in subdomains of intermediate size tends to the exact value 
on increasingly coarser levels. The accuracy of averages (including the particle number 
fluctuation) depends only on the CP tables adequacy, which is achieved by the correct 
selection of the variable set. 

The advantage of the multilevel Monte Carlo method consists in fast convergence of 
measured mean values due to the selfconsistent equilibration on different levels modes of 
wave lengths comparable with the simulation domain. The computational work on each 
level is proportional to the number of gridpoints and is independent of the particle number 
associated with the gridpoint. It leads to the high speed of the method as compared with 
conventional algorithms. 

The equilibrium on fine levels (and in particular the frequency of the appearance of 
a given particle number in the simulation domain) is defined by the canonical ensemble 
configurations on coarsest level. The particle number in fine level simulation domains is 
variable and its distribution imitates the result of the grand canonical ensemble simula- 
tion 120]. The multilevel Monte Carlo method is not concerned with the equilibrium in 
accordance with the value of a chemical potential. It opens new way for the development 
of the one-phase approach to the phase-transition problem fIT\ . 
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